Benchmarks

Numerical Improvement

I created a test where I initialize a boundary at some specified resolution and nfp and run both the nfp optimized code and with the same boundary but with nfp = 1 after I mirror the boundary such that it goes through the dense solve of the code. There is clearly a great speedup/memory reduction that increases with nfp, showing I am taking advantage of the symmetry. Note that these do have AdaptiveArrayPooling integrated but were tested with MAYBE_POOLING[] = false, a cool testing feature Min-Gu showed me that allows you to turn it off for selective functions to see what the memory usage looks like.

small  (N=1152)  N=3456 (24×48/period × nfp=3)  M=1152  modes=15  residue-classes=3
------------------------------------------------------------------------------
  full (dense N×N)        time= 0.769 s  @alloc=  193 MiB
  reduced (block-circ.)   time= 0.306 s  @alloc=  109 MiB
  speedup: 2.51×   max rel err: 1.51e-15

medium (N=2304)  N=6912 (48×48/period × nfp=3)  M=2304  modes=42  residue-classes=3
------------------------------------------------------------------------------
  full (dense N×N)        time= 6.549 s  @alloc=  769 MiB
  reduced (block-circ.)   time= 2.246 s  @alloc=  453 MiB
  speedup: 2.92×   max rel err: 1.59e-15

large  (N=4608)  N=9216 (48×48/period × nfp=4)  M=2304  modes=42  residue-classes=4
------------------------------------------------------------------------------
  full (dense N×N)        time=14.378 s   @alloc= 1351 MiB
  reduced (block-circ.)   time= 3.068 s  @alloc=  521 MiB
  speedup: 4.69×  max rel err: 3.80e-15

3D GPEC Integration

For an integration test with 3D GPEC + n_stride capabilities, I ran the full workflow for Junyoung's sample NFP = 3 equilibrium with:

n_low = 1, n_low = 6, n_stride = 1 -> n = [1, 2, 3, 4, 5, 6] (previously working case)
[RESULT] dW_total  = 6.608533e-01
[RESULT] dW_plasma = 2.212641e-01
[RESULT] dW_vacuum = 4.395892e-01

and

n_low = 1, n_low = 6, n_stride = 3 -> n = [1, 4]
[RESULT] dW_total  = 6.436734e+01
[RESULT] dW_plasma = 1.653465e+00
[RESULT] dW_vacuum = 6.271387e+01

and

n_low = 2, n_low = 5, n_stride = 3 -> n = [2, 5]
[RESULT] dW_total  = 6.608515e-01
[RESULT] dW_plasma = 2.212649e-01
[RESULT] dW_vacuum = 4.395866e-01

This result indicates that the most unstable mode lives in the n = [2, 5] decoupled family, and we can recapture it to within numerical error while only running with these modes. When running with the other n = [1, 4] family, we capture a more stable mode instead

@claude review

Review

This is a large PR (710+/1582-, 18 files) that bundles two things: (1) the NFP block-circulant fast path for 3D vacuum response, and (2) a broad refactor of the Fourier-transform plumbing from grouped real/imag (cos/sin) matrices to a single complex basis, which ripples through ForceFreeStates, PerturbedEquilibrium, and ForcingTerms. The author already flagged in the description that (2) probably should have been its own PR — agreed, but understandable given the dependency.

Bugs / correctness

• Duplicate dead code in _compute_vacuum_response_3d! (src/Vacuum/Vacuum.jl): the destructuring and zeroing lines are repeated back-to-back:

  (; mtheta, nzeta, nfp, m_modes, n_modes) = inputs
  fill!(vac_data.wv, 0)
  
  (; mtheta, nzeta, nfp, m_modes, n_modes) = inputs
  fill!(vac_data.wv, 0)

  Harmless (idempotent), but it's an obvious copy-paste leftover that should be collapsed to one copy before merge.

• Orphaned structs: KernelParams2D / KernelParams3D in src/Vacuum/DataTypes.jl are no longer constructed or consumed anywhere — the kernel! dispatch wrappers that used them were removed from Kernel2D.jl/Kernel3D.jl in favor of calling compute_2D_kernel_matrices!/compute_3D_kernel_matrices! directly with raw arguments. Worth deleting the now-dead structs (and their docstrings) rather than leaving them to rot.

Physics / numerics

• The block-circulant reduction (_compute_vacuum_response_3d!, nfp > 1 branch) looks structurally correct: the per-residue-class accumulation D̂ₖ = Σ_d D_d·ω^{kd} matches the docstring, the block indexing into grad_green/green_temp lines up with the field-period-tiled source columns, and the final mode-space projection via E * (G_plasma * E') is consistent with how the 2D path projects via ft.basis. Good that this is backed by a real test (compute_vacuum_response 3D field-periodic reduction) that checks the reduced block-circulant path reproduces the independently computed dense full-torus result (nfp=1 path) to rtol=1e-6, plus a residue-class decoupling check and Hermitian-symmetry check. That's a meaningful self-consistency check given there's no Fortran reference for this feature.
• grri/grre are correctly zeroed and explicitly documented as unsupported in the 3D path — fine as a documented gap, but worth double-checking no caller in PerturbedEquilibrium could silently pick up zeroed 3D Green's functions once 3D singular-coupling work starts landing.
• The real/imag → complex Fourier basis rework (FourierTransforms.jl, CoilFourier.jl, ResponseMatrices.jl, SingularCoupling.jl) is consistent end-to-end as far as I traced it — shapes, transpose conventions, and the dropped factor of 2 on grri/grre (now single complex array instead of grouped [2N, 2*mpert]) all line up with the new mpert × mtheta basis convention. FourierTransforms.jl shrinks from ~700 to ~130 lines by dropping the redundant in-place/grouped-real-imag API, which is a real simplification, not just churn.

Process

• Regression harness: per CLAUDE.md, this is mandatory before merging any substantive change, and this PR is about as substantive as it gets — it touches CoordinateInvariant.jl, ForceFreeStatesStructs.jl, Free.jl, ResponseMatrices.jl, and SingularCoupling.jl, all of which feed the existing 2D regression-tracked pipeline (diiid_n1, solovev_n1, solovev_multi_n), even though the headline feature is 3D-only. I don't see a regression report in the PR description — please run regress --cases diiid_n1,solovev_n1,solovev_multi_n --refs develop,local and post the table before merging. Given there's no existing 3D regression case, it's also worth a follow-up to add one once grri/grre support lands for the 3D path.

Minor

• Deleting the stale benchmarks/benchmark_fourier_transforms.jl is reasonable cleanup given the API it benchmarked no longer exists.
• A fair amount of the diff is docstring reflow (bullet indentation/blank-line changes) unrelated to the functional change — not wrong, but it inflates the diff and makes the substantive parts harder to spot in review. If this was an auto-formatting pass, consider calling that out explicitly in the PR description next time.

Regression results for solovev_n1, solovev_multi_n, diiid_n1
Summary: 21 unchanged
Summary: 15 unchanged
Summary: 46 unchanged

@jy021125 wanted to make sure you are tagged on this one - this should enable the n_stride = nfp implementation we discussed over email + some runtime improvements to speed and memory. I have pushed the corresponding changes needed in the python to the open PR on your github, so if you want to test this, it should work with the python on the vacuum_interface branch and the julia on this PR's branch.

Let me know if you find any bugs testing this out!